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Summary 

In this work we deal with correlated failure time (age at onset) data arising from population- 
based case-control studies, where case and control probands are selected by population- 
based sampling and an array of risk factor measures is collected for both cases and controls 
and their relatives. Parameters of interest are effects of risk factors on the hazard func- 
tion of failure times and within-family dependencies of failure times after adjusting for 
the risk factors. Due to the retrospective nature of sampling, a large sample theory for 
existing methods has not been established. We develop a novel estimation techniques for 
estimating these parameters under a general semiparametric shared frailty model. We 
also present a simple, easily computed, and non-iterative nonparametric estimator for the 
cumulative baseline hazard function. A rigorous large sample theory for the proposed 
estimators of these parameters is given along with simulations and a real data example 
illustrate the utility of the proposed method. 

Keywords 

Case-control study; Correlated failure times; Family study; Frailty model; Multivariate 
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1 Introduction 



Clustered failure times arise often in medical and epidemiologic studies. Examples include 
disease onset times of twins (in terms of age), multiple recurrences times of infections on an 
individual, or time to blindness for both eyes within an individual. A typical case-control 
family study includes a random sample of independent diseased individuals (cases) and 
non-diseased individuals (controls), along with their family members. An array of genetic 
and environmental risk factor measures is collected on these individuals. Integration of 
genetic and environmental data is a central problem of modern observational epidemiology 
(Hopper et al., 1994; Zhao et al., 1997; Malone et al., 1998; Malone et al. 2000; Becher et 
al., 2003). Case-control family studies are powerful because they provide an efficient way 
to assess the effect of risk factors on the occurrence of a rare disease, and furthermore 
allow researchers to dissect genetic and environmental contributions to the disease based 
on the familial aggregation pattern of disease clusters. Hopper (2003), in a commentary 
article, suggested that such study designs may be the future of epidemiology, not just 
genetic epidemiology. Hence the need for statistical methods that can fully utihze such 
data is acute. 

In this work we focus on population-based case-control family studies, where a number 
of case and control probands are randomly sampled from a well-defined population. The 
probands are the index subjects because of whom the families are ascertained. Here we 
use the term proband in a broad sense to refer to both cases and controls, in contrast 
with the traditional usage in which proband refers only to cases. 

Relative to classical case-control methods, analysis of such studies is complicated in 
several ways: (1) Comparisons are no longer solely between subjects with and without the 
disease under study, but rather between collections of the case probands and their relatives 
and the control probands and their relatives, each collection typically including many 
subjects both with and without the studied disease. (2) Data are clustered within families, 
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and hence reflect intra-familial correlation due to unmeasured genetic and environmental 
factors. 

Our work is motivated by a recent breast cancer study conducted at the Fred Hutchin- 
son Cancer Research Center (Malone et al., 1998; Malone et al., 2000). In this study, the 
cases were incident breast cancer cases ascertained from the Surveillance, Epidemiology, 
and End Results (SEER) registry, which is a set of geographically defined, population- 
based cancer registries in the United States. The controls were selected by random digit 
diahng, and were matched with cases based on age at diagnosis and county of residence. 
Female relatives of case and control probands were identified, and the risk factor and 
outcome information was subsequently collected on these relatives. The primary goals of 
the study are (a) to determine the degree of the strength of the dependency of ages at 
diagnosis of breast cancer between probands and their relatives; (b) assess the effects of 
covariates on breast cancer risk. 

Two modeling approaches, marginal and conditional, are typically used for accounting 
for the correlation within a cluster. In the conditional model, the correlation is explicitly 
induced by a cluster-specific random effect, with the outcomes of the cluster members 
being conditionally independent given the random effect. The random effects model for 
failure time outcome is generally known as frailty model, in that the random effect or 
frailty is assumed to act multiplicatively on the baseline hazard rate of failure. Many 
frailty models have been considered, including gamma (Gill, 1985, 1989; Nielsen et al., 
1992; Klein 1992, among others), positive stable (Hougaard, 1986; Fine et al., 2003), 
inverse Gaussian, compound Poisson (Aalen, 1992) and log-normal (McGilchrist, 1993; 
Ripatti and Palmgren, 2000; Vaida and Xu, 2000, among others). Hougaard (2000) 
presented a comprehensive review of the properties of various frailty distributions. Under a 
frailty model, the regression coefficients are cluster-specific log-hazard ratios. By contrast, 
in the marginal model the correlation is modelled through a multivariate distribution, such 
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as a copula function (Genest and MacKay, 1986; Marshall and Olkin, 1988; Shih and Louis, 
1995) with a specified model for the marginal hazard functions. The regression coefficients 
in the marginal model represent the log-hazard ratios at the population level regardless 
of which cluster an individual comes from. The effect therefore is "population-averaged." 
Zeger et al. (1988) provided a comprehensive comparison of the conditional and marginal 
modelling approaches. 

Methods have been developed for the age at disease onset data from case-control family 
studies under both modelling approaches. Shih and Chatterjee (2002) proposed a semi- 
parametric quasi-partial-likelihood approach for estimating the regression coefficients in a 
bivariate copula model. Their cumulative hazard estimator requires an iterative solution, 
and thus the properties of their estimators could only be investigated so far by a simulation 
study. Moreover, in the presence of multiple relatives for each proband, the relatives were 
treated as if they were independent of each other, which may lead to loss of efficiency 
in the baseline hazard function estimator. In contrast, Hsu et al. (2004) presented a 
quasi-EM algorithm method for the popular gamma frailty model. In the random effects 
model of Hsu et al., the regression coefficients express the effect on a subject's disease risk 
due to being exposed relative to the same subject's level of risk when unexposed. The 
baseline hazard function estimator naturally accommodates multiple relatives in a family 
(Hsu and Gorfine, 2006). However, the properties of the proposed estimators were also 
studied only by simulation. The method of Shih and Chatterjee (2002) can be adapted to 
the family-specific frailty setting (Oakes, 1989), but with the same limitation as for the 
marginal model: the lack of large sample theory. 

In this work, we develop a new estimation technique for the general semiparametric 
shared frailty model, where the parameters of interest are the regression coefficients and 
the frailty parameters. Our general family-specific frailty model is for any frailty distribu- 
tion that has finite moments. The estimation procedure for the baseline hazard function 
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leads to an estimator whose asymptotic properties can be derived and expressed in a 
tractable manner. 

Section 2 presents our model, and Section 3 describes our estimation procedure. Sec- 
tion 4 gives the consistency and asymptotic normality results for the estimators. In 
Section 5, we describe an extension of our method for the case where the proband obser- 
vation times are subject to a certain restriction that can arise in some studies. Section 
6 presents simulation results. In Section 7 we illustrate our method with a case-control 
family study of early onset breast cancer. Section 8 provides a short discussion. The 
Appendix provides the details of the asymptotic theory. 

2 Notation and model formulation 

We consider a matched case-control family study where one case proband is age-matched 
with one control proband, and an array of risk factors is measured on the case and 
control probands and their relatives. Each matched set contains one case family and one 
control family, and there are n i.i.d. matched sets. Let T?- and Cjj denote the age of 
disease onset and age at censoring, respectively, for individual j of family i, i — 1, . . . , 2n, 
j — 0,1, ... , rrii, where j — corresponds to the proband. Following Parner (1998, p. 187), 
we regard mi as a random variable over {1, . . . , m) for some m, and build up the remainder 
of the model conditional on m^. Define 5ij — I{T^j < Cij) to be the failure indicator and 
Tij — min(T°-, Cij) to be the observed follow-up time for individual ij. We assume that a 
p- vector of covariates is observed on all subjects, and let Z^j denote the value of the (time- 
independent) covariate vector for individual ij. In addition, we associate with family i 
an unobservable family- level covariate cui, the "frailty", which induces dependence among 
family members. The conditional hazard function for proband i, given the family frailty 
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uji^ is assumed to take the form 

Aio(i|Zjo,Ci;i) =u;jAo(t)exp(/3^Zjo) i = l,...,2n. (1) 

The conditional hazard function for relative ij, j — 1, . . . , mj, given the family frailty uji 
and proband z's data, is assumed to take the form 

Aij(t|Tjo,5jo,Zio,Zij,a;j) =u;jAo(t)exp(^^Zjj) i = l,...,2n; j = l,...,mj. (2) 

Here /3 is a p-vector of unknown regression coefficients, and Aq is a conditional baseline 
hazard of unspecified form. The above model implies that the proband and the relatives 
have a common conditional baseline hazard function Aq, and that all the dependence 
between the proband and the relatives in a given family is due to the frailty factor uii. 
The random variable uji is assumed to have a density f{uj) = f{w;9), where 9 is an 
unknown parameter. For simplicity, we assume that ^ is a scalar, though the vector case 
could be developed in a similar manner. 

We put 7 = (/3^, 9)^ , and let 7° = (/3°^, 9°)^ denote the true value of 7. The objective 
is to estimate 7 and Ao{t) — Xo{u)du. Let A.Q{t) denote the true value of Aq. Further, 
let Siu — {Sii, . . . , 5imi), Tiji — (Til, . . . , Tj„.), and Zj^ = (Zji, . . . , Zj^.). 

We make the following basic assumptions. 

1. Zij is bounded. 

2. The parameter 7 lies in a compact subset Q of IR^"*"^ containing an open neighbor- 
hood of 7°. 

3. Conditional on {Zij}^^^ and Ui, the censoring times are independent and noninfor- 
mative for cui and Aq). In addition, the frailty cui is independent of {Zjj}^^. 

4. The effect of the covariates on age at onset is subject-specific, i.e. Pr(Tjj , Sij \ Zjo, Zjj^, cui) 
= Pr{Tij, 6ij\Zij,u;i). This implies Pr{Tij, 5jj|Zjo, Zj^j) = Pr{Tij, Sij\Zij). 
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The first two of tliese assumptions imply tliat tliere exists a positive constant v sucli 
tfiat 

v'^ < exp(/3^Zij) < u. (3) 

A number of additional technical assumptions are listed in the appendix. 
The likelihood function for the data can be written as 

2n 

iR, "^iR, Zjo|Tjo, Sio) 

1=1 
2n 

= Zjo,Tio,5io) X /(Zji?|Zio) X /(Zjo|Tjo,(5io)- (4) 



1=1 



Since /(Zjij|Zjo) does not depend on the parameters of interest (/3,Ao,6'), this term will 
be ignored. In the following subsections we consider the other two terms in (jl]). 

2.1 The likelihood for the proband data 

For the likelihood function of the proband data, fl^" /(Zjo|Tjo, 610), we use a retrospective 
likelihood for the standard case-control study (Prentice and Breslow, 1978). We express 
this likelihood in terms of the marginal survival function Sio{t) = Pr{Tio > t|Zjo) = 
J Pr{TiQ > t\ZiQ,uj)f{u!)duj. In our setting we have n one-to-one matched sets. Based on 
the marginal survivor function, the marginal hazard function can be written as 

Ko{t) = Ao(t)exp(/3 Zio) — — , 

where 

/ifci(t;7, Ao) = j uj^ exp{-ujHio{t)}f{uj)duj /c = 0,l,2 

and HiQ(t) = Ao(t) exp(/3"^Zjo). We arrange the notation so that the first n families are the 
case families and the rth case family, r = 1, . . . , n, is matched with the {n + r)th control 
family. The likelihood for the proband data is then replaced by the following conditional 
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likelihood: 



^(1) ^ TT exp(/3^Z^o)^io,.(T,.o; 7, Ap) 

fJl exp(/3^Zro)ClOr(rrO; 7, ^0) + exp(/3'^Z(„+r)o)6o(n+r)(7'(n+r)o; 7, ^o) ' 



where 



Let 



^fefe/j(i;7, Ao) = 7- 7-^, A;,A; =0,1,2. 



eS,(t;7,Ao) = ^eioi(t;7,Ao) = ^ioW^io/ 7, Ao) - e20i(t; 7, Ao)} , 



and 



d 

Cfcfc'i(^;7,Ao) = ^CfcA:'i(t;7,Ao) 

^ /Xfc'i(t; 7, ^o)l^kii't'^ 7, Aq) - 7, Ao)/Xfc^(t; 7, Ap) 

/^fc'i(^;7,Ao) 

for /c, /c' = 0, 1, 2 and I — 1, ... ,p. Then the score function for Pi, I — 1, ... ,p, is given by 



r=l 



?10r (7;o;7,Ao) 

_ cxp(/3^Z^p)[Z^o;^io^(T,.o; 7, Aq) + Cfo^(T,.p; 7, Ap)] 

exp(/3^Z,.o)^iPr.(T^o; 7, Ap) + cxp(/3^Z(„+^)o)^io(„+r)(T(„+r)p; 7, Aq) 

exp(/3^Z(„+r)o)[%t+r)0«'ClO(n+r)(T(„+,.)p; 7, Ap) + 'Cl0(„+r)(^(n+r)0; 7> Aq)] 

exp(^^Zro){io^(Tro; 7, Ap) + exp(^^Z(„+^)p){ip(„+r)(T(„+^)o; 7, Ap) 
and the score function for 9 is given by 



l?10rUrO,7)-'»-0 



^ l6or(T^o;7) Ap) 

exp(/3^Z^,o)^fo,.(T^o; 7, Ap) + exp(/3^Z(„+,.)p)^^Q(^^^^)(T(„+^)o; 7, Ap) 
exp(^^Z^p)^ip^(T^p; 7, Ap) + exp(;3^Z(„+^)p)^ip(„+^)(T(„+^)o; 7, Ap) 

Under the gamma frailty model, we have 

/XH(t;7,Ao) _ rflu- , 
//pi(t, 7, AqJ 
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and so the likelihood function corresponds to that presented in Hsu et al. (2004) in the 
case of one-to-one matching. Extension to matching of multiple cases or multiple controls 
are straightforward, see e.g. Breslow and Day (1980). 



2.2 The likelihood for the data from the relatives 

Let Nij{t) = Si,I{Ti, < t), J = l,...,m„ iV,.(t) = ETliNij{t), Hi,{t) = Ao(T,,- A t) 
exp(/3"'"Zjj), j = 1, . . . , rrii, and Hi {t) = Yljli Hij{t), and let r be the maximum follow- 
up time. The likelihood of the data from the relatives then can be written as 

2n ^ m-i 

ri /n-iln-i [- r r \-i(5^o^/n-iln-i r r * \ p / Im T r-w \ ^ 

w 



JJ^ {\ij{Tij\TiQ, (5io, Zjo, io)} Sij{Tij\Tio, 5jo, Zjo, uj)f{w\Tio, 5^0, Zio)(i' 

3=1 

2n nii 2n „ 

n n exp(/3^Z,,)}''^ J] / exp{-a;i/,.(r)}/(cu|r,o, <5,o, Z,o)rfcu. 



i=l j=l 1=1 

Here, by a Bayes theorem argument 



<5iO 



exp(-a;Ao(T,o)e^"z»)/(a;) 



f (c(j|T,n, (5,ni Z,n) = ^- . (6) 



The log-likelihood is given by 

2n rrii 2n 



^^'^ = $^X]'^^^-^°g{^o(T.,)exp(/3'^Z,,)}+J]log j / a;^> Mexp{-cuff,,(^)}/(^l^.o,5io,Zio)rfcu 

1=1 j=l i=l ^"^ 



The scores for (/3i, . . . , /3p) are given by 
r(2) 



/^^-Mexp{-^iJ,,(^)}4/(^l^io,5io,Z,o)rfa; 



' °^ " trfe tr /^^'•(-)exp{-a;ff,.(r)}/(c.|T,o,5.o,Z.o)rfc. 

^ / exp{-^g,.(^)}/(^|7^^o, ^^0, ^io)duJ A 

4- /.;^«.Mexp{-^i7,.(^)}/(^|7;o,5^o,Z.o)ci^ ft ' 

2 — 1 7 1 



for / = 1, . . . ,p. The score for 9 is given by 



..(2) . A N _ / uj^'-^^'> exp{-ujHiXT)}^f{uj\Tio, 6io, Zio)du; 
p+il7, /a;iv..(r)exp{-u;i7,,(^)}/(^l^.o,^.o,Z,o)cia; 
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3 The proposed approach 



We focus first on estimating the baseline cumulative hazard function Ao(t). Let Yij{t) = 
I(Tij — ^)) let J-'t denote the cr-algebra generated by (Tjo, (5jo, Zjo) plus the entire 
observed history of the relatives up to time t: 

Tt = a{Tio,5io,Zio,Nij{u),Yij{u),Zij;i = 1, . . . , 2ra; j = 1, . . . , m^; < m < t). 

Then, as discussed by Gill (1992) and Parner (1998), the stochastic intensity process for 
Nij(t), i = 1, . . . , 2n, j = 1, . . . , rrii, with respect to J-'t is given by 

Ao(t) exp{f3^Zij)Yij{t)^i{t-, 7, Aq), (7) 

where, using ([H]), 

V'i(t,7,Ao) = E[uJi\J='t] 

J ^^-.W+i exp{-ujHat))f{uj\T,o, 5^o, Z,o)duj 
J exp{-uHdt))f{oo\T,o, 6,0, Z,o)du 

J ^A., exp(-cu{iJ,.(^) + H,o{T,o)})f{uj)duj 

Define (for < r < m and h > 0) 

Jw^+'e-'^^f{w)dw 
^ ^' ^ = r r -hwf( w • ^ 

Some salient properties of ip*{r,h) are noted in Sec. 9.2. With this definition, we have 

The key to obtaining parameter estimators for a semiparametric survival model is 
an estimator of the nonparametric baseline hazard function. For our model, a Breslow- 
type estimator with a jump at each observed failure time among the relatives can be 
formulated in a natural way (Shih and Chatterjee, 2002). However, the hazard function 
for the relatives at time t depends on their respective proband's observation time TjQ. For 
example, under the gamma frailty model with expectation 1 and variance 6, tpi{t, 7, Aq) = 
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{6~^ + Ni^t) + Sio}{9~^ + Hi it) + iJjo(^io)} Often the relevant proband's observation 
time is greater than t, so that the standard Breslow formula for the baseline hazard 
estimator at time t involves values of Aq for times beyond time t. An iterative procedure is 
thus required to obtain the estimator. In addition, because of this estimator's complicated 
structure, its asymptotic properties have not been established. 

We propose to estimate the baseline hazard function using a non-iterative two-stage 
procedure. The first-stage estimator is a weighted Breslow-type estimator, where the 
weight at time t for family i is equal to 1 if the observation time Tjo of the family i 
proband is less that t, and equal to otherwise. The second-stage estimator is the 
standard Breslow-type estimator that uses all the relatives' failure times, plugging in the 
first-stage estimator where necessary. 

More specifically, the estimators are defined as follows. Let r^,, g = 1, . . . , G, denote 
the observed failure times of the relatives and assume that dg failures were observed at 
time Tg. In theory, since we are dealing with continuous survival distributions, dg = 1 for 
all g, but we write the formula for the estimator in a form that allows for a modest level of 
ties in the survival times. Let Amax be some known (possibly large) upper bound for Ag(t). 
Define il){r^h) = ip*{r,h A hmax), with hmax = fn^Amax, where z/ is as in ([3]). Further, 
define 7,A) = ^/'(A^j(t), 7, A)). The first-stage estimator is then defined as a 
step function whose g-th jump is given by 

Ef:,mo<rg)j:-l,dN,,irg) 

' Eti mo < rgmr,-in, Ao) E^i Y.^r,) exp(/3^Z,,) ' 

In a similar way, the second-stage estimator is defined as a step function whose g-i\i jump 

is given by 

AAo(rg) = — — . , (10) 

Ei=i A{.rg-u 7) Er=i ^^P(/^ ^^i) 

where ipiit^l) is defined analogously to ipiit, ^, Aq), with Ao(Tjo) replaced by Ao(Tjo) if 
Tio > t and by Ao(Tjo) otherwise. It is clear that no iterative optimization process is 
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required here and the large-sample properties of Ao(i) will be determined by those of 
Ao(t). 

We note that there is no guarantee that Ao(i,7) as defined above will be bounded by 
A-max-i but this does not matter: if desired, we can replace the estimator by min{Ao(t, 7), A^a^;} 
without affecting the asymptotics. 

For estimating (/3, 0) we use a pseudo-likelihood approach: in the score functions based 
on L^^) and L^^^ , we replace the unknown Aq by Aq. Thus, the score function corresponding 
to A (for / = 1, . . . ,p) is given by [7^7, Ao) = n'^ |[//'^(7, Aq) + t/f ^(7, Aq)} , and the 
estimating function for 9 is given by ^'^+1(7, Aq) — |t/p^\(7, Aq) -|- C/p^\(7, Ao)| . To 
summarize, our proposed estimation procedure is as follows: 

1. Provide an initial value for 7. 

2. For the given values of 7, estimate Aq using (9) and (10). 

3. For the given value of Aq, estimate 7. 

4. Repeat Steps 2 and 3 until convergence is reached with respect to Aq and 7. 

4 Asymptotic properties 

We show that 7 is a consistent estimator of 7° and that \/n(7 — 7°) is asymptotically 
mean-zero multivariate normal. In this section, we present a broad outline sketch of the 
argument. The Appendix provides the details of the proofs, including a detailed hst of 
the technical conditions assumed. The arguments are patterned after those of Gorfine et 
al. (2006) and Zucker et al. (2006), but with considerable expansion. 
Consistency is shown through the following steps. 

Claim Al. Ao(t, 7) converges in probability to some function AQ(t, 7) uniformly in t and 
7. The function A^(i,7) satisfies A*(t,7°) = A°(i). 
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Claim A2. Ao(t, 7) converges in probability to some function Ao{t, 7) uniformly in t and 
7. The function Ao(i,7) satisfies Ao(i,7°) = Ao(^)- 

Claim B. U(7, Ao(-, 7)) converges in probability uniformly in t and 7 to a limit u(7, Ao(-, 7)). 

Claim C. There exists a unique consistent (in pr.) root to U(7, Ao(-,7)) = 0. 

The proofs of Claims Al, A2, and B involve empirical process and function-space com- 
pactness arguments, while Claim C is shown using Foutz's (1977) theorem on consistency 
of maximum likelihood type estimators. 

The proof of asymptotic normality is based on the following decomposition: 

0= U(7,Ao(-,7)) 

= U(7°, A°) + [U(7°, Ao(-, 7°)) - U(7°, A°)] 
+ [U(7,Ao(-,7))-U(7°,Ao(-,7°))]- 

In the Appendix we analyze each of the above three terms and prove that -\/ri('7 ~ 7°) 
has an asymptotic mean-zero multivariate normal distribution. Although it is possible to 
develop a consistent closed-form sandwich estimator for the asymptotic covariance matrix 
of this distribution, we do not present this estimator because it is too complicated to be 
practically useful. Instead, as discussed in Section 6, we recommend bootstrap standard- 
error estimates. 

5 Extension to restricted sampling of probands 

A key assumption in our procedure for estimating Aq is that the support of the proband 
observation times and that of relatives' observation times have the same lower limit, which 
is designated (without loss of generality) as time zero. In some applications, however, the 
probands' observed times are restricted to some range [sq, Si] with sq > 0. For example. 
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Malone et al. (2006) present a multi-center case-control breast cancer study where ages 
of cases and controls are restricted between ages 35-64. In a design of this form, where 
the probands' observed times are left-restricted by sq and the relatives' failure times are 
unrestricted, Aq will be underestimated by our two-stage procedure. But this bias can be 
easily corrected by first estimating Ao(so). 

We present here the resulting three-stage estimator for the left-restricted design. Let 
AAo{rg, Ao(so)} and AAofr^, Ao(so)} be defined analogously to AAo(rg) and AAo(rg) 
with Ao(Tjo) = Ao(so) + J^Tg&iso Tm] ^o(tc,)- The estimator Ao(so) is defined to be the root 
of 

J2 AAo{r„Ao(so)}-Ao(so) = 0. (11) 

rgG[0,so] 

The root can be found by simple univariate Newton-Raphson iteration. This completes 
the first stage. The second stage involves calculating AAqIt^, Ao(so)}, g = 1, . . . ,G, using 
the formula (Q. In the third stage, we use the results of the second stage and the formula 
( fTOl) to calculate the the final estimate AAo^Tg), g = 1, . . . ,G. In applying (fTOll . we replace 
Ao(Tio) by Ao{Tio, Ao(so)} if T^o > Tg and by Ao(T;,o) otherwise. 

In Section 6 below, we present simulation results for this estimator. In theory, the 
asymptotic properties of the three-stage procedure could be worked out via an extension 
of the arguments for the two-stage procedure, but the algebra becomes very complicated. 
We hope to develop asymptotic theory for the left-restricted design in future work. 

6 Simulation results - gamma frailty 

We have performed a simulation study to evaluate the finite sample performance of the 
proposed method and compare it with existing methods. One of the most extensively used 
frailty models is the model with gamma-distributed frailty. Under this model, 6 quantifies 
the heterogeneity of risk among families. The larger the value of 6 is, the stronger the 
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dependence among family members. In addition, the gamma frailty model can be re- 
expressed in terms of the Clayton-Oakes copula-type model (Clayton, 1978; Oakes, 1989) 
and the cross-ratio, introduced by Oakes (1989) as a local measure of association between 
survival times, is constant on the support of failure time region and equals 1 + 9. The 
gamma frailty model is also convenient mathematically, because it admits a closed-form 
representation of the marginal survival distributions. These features make the gamma 
frailty model very popular. Hence we conducted our simulation study under the gamma 
frailty model, using, as is customary, the gamma distribution with expectation 1 and 
variance 9. 

Simulation results are based on 500 control probands matched to 500 case probands, 
with one relative sampled for each proband. We considered a single U[0, 1] distributed 
covariate with (3 = ln(2), Ao(t) = t, 9 = 2, and a t/[0, 1] censoring variable, yielding a 
censoring rate among the relatives of approximately 60%. In Table 1 we compare the 
following three estimates: the proposed estimate with the two-stage procedure for Aq, the 
estimate of Hsu et al. (2004), and a modified version of Shih and Chatterjee's (2002) 
estimate, with their method adapted to the gamma frailty model. Results are based on 
500 simulated data sets. The efficiency difference between our two-stage estimator and 
that of Shih and Chatterjee is very small. 

For our estimators, in addition to the above-mentioned simulation setting, we also 
considered /3 — 0, 9 — 3 and a censoring distribution of C/[0, 4] with a C/[0, 4] distributed 
covariate, or a censoring distribution of U[0, 0.1] with U[0, 1] distributed covariate, yield- 
ing censoring rates of approximately 30% or 90%, respectively. To construct confidence 
intervals, we use a bootstrap approach. In the setting of censored survival data, the 
usual nonparametric bootstrap is problematic because it leads to a substantial proportion 
of tied survival times. Hence we used the weighted bootstrap approach of Kosorok et 
al. (2004) instead. For the weighted bootstrap, a sample of 2n independent and iden- 



16 



tically distributed weights from the unit exponential distribution was generated for each 
bootstrap sample. Let ^i, . . . , ^2n be the standardized weights after dividing each weight 
by the average weight. Then, in the estimating functions, for any given function h the 
empirical mean Ylf^i h{Ti, 5i, Zj) is replaced by its corresponding weighted empirical 
mean n~ Ei=i ^MT^, S„ Z,). Kosorok et al. (2004) proved that this weighted bootstrap 
procedure gives valid inference for all parameters under right- censored univariate failure 
times. 

Results based on the two-stage procedure for Aq are presented in Tables 2-4 for various 
levels of censoring. We present the mean, the empirical standard error, and the coverage 
rate of the 95% weighted bootstrap confidence interval. The results are based on 50 
bootstrap samples for each of the 2000 simulated data sets of each configuration. Our 
estimates perform well in terms of bias and coverage probability. 

For studying the case of left-restricted data, we considered a similar configuration as 
of Table 1, but now the probands observation times are restricted to be > 0.1. In Table 
5 the results of our three-stage estimator are presented along with the estimators of Hsu 
et al. (2004) and of Shih and Chatterjee (2002). It is seen that estimating Ao(so) yields 
small efficiency loss in Aq, in compare to the other two methods. 

7 Example 

We apply our method to the breast cancer study mentioned in the introduction. Various 
risk factors were measured on probands and their relatives. For illustrative purposes we 
consider age at first full-term pregnancy with the relatives of the probands being the 
mothers. The following analysis is based on 437 breast cancer case probands matched 
with 437 control probands and a total of 874 mothers. The number of mothers who 
had breast cancer was 70 among the case families and 35 among the control families. 
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The number of women whose first five birth occurred before age 20 was 142 among the 
probands and 181 among the mothers. In the following analysis, the gamma frailty model 
is used with expectation 1 and variance 9. Three estimation procedures are considered: 
our proposed method, the Hsu et al. (2004) method, and the modified Shih-Chatterjee 
(2002) method. For our proposed method, the two-stage procedure for Aq is used since 
the age range of the mothers with breast cancer was 20-76 and of the age range of the 
probands was 22-44. Table 6 presents the regression coefficient parameter estimate 
the dependency parameter estimate, 6, and Aq at ages 40, 50, 60 and 70 years old, along 
with their respective bootstrap standard errors. The proposed approach and that of 
Shih and Chatterjee yielded similar dependency estimates with the proposed approach 
being moderately more efficient. Hsu et al.'s approach gave a slightly lower dependence 
estimate. The regression coefficient estimates of Hsu et al. and that of Shih and Chatterjee 
are similar, with the latter being slightly more efficient. The proposed approach yielded 
a slightly lower covariate effect. The cumulative baseline hazard estimates are similar 
under the three estimation techniques. The results, based on the three methods, imply 
that women who had their first full-term pregnancy before age 20 have a reduced risk of 
developing breast cancer, supporting the observation of breast cancer risk reduced by early 
first full-term pregnancy (e.g. Coditz et al., 1996; among others). The estimates of the 
dependency parameter imply that after adjusting for the first full-term pregnancy, there 
remains a significant dependency between the ages of onset for mothers and daughters 
with cross ratio (1 + 6*) close to 2. 

8 Discussion 

In this work we have presented a new estimator for case-control family study survival 
data under a frailty model, allowing an arbitrary frailty distribution with finite moments. 
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Rigorous large sample theory has been provided. Simulation results under the popular 
gamma frailty model indicate that the proposed procedure provides estimates with min- 
imal bias and confidence intervals with the appropriate coverage rate. Moreover, our 
estimators were seen to be essentially identical in efficiency to estimators based on the 
more complex approach of of Shih and Chatterjee (2002). 

Rigorous large sample theory has been provided for unrestricted sampling of probands. 
For restricted sampling, the asymptotic theory could be worked out largely following the 
arguments for the two-stage estimator but the algebra becomes very complicated. It is 
beyond the scope of the current paper and will be presented in a future communication. 

9 Appendix: Asymptotic theory 

This appendix presents the technical conditions we assume for the asymptotic results and 
the proofs of these results. The development is patterned after Zucker (2005) and Zucker 
et al. (2006), but considerable extension of the arguments is required. In the presentation 
below, we focus on the added arguments needed for the present setting, and refer back to 
Zucker (2005) and Zucker et al. (2006) for the other segments of the development. 

9.1 Assumptions and background 

In deriving the asymptotic properties of 7, we make a number of assumptions. Several 
of these assumptions have already been presented in the main text. Below we list the 
additional assumptions. 

1. There is a finite maximum follow-up time r > 0, with ^[Yl^ii^iji'T)] — V* > ^ for 
all i. 

2. The frailty random variable cui has finite moments up to order {■m + 2). 



19 



3. There exist 6 > and C > such that 

hmw-(*-^)/(uO = C. 



4. The basehne hazard function X^lt) is bounded over [0, r] by some fixed (but not 
necessarily known) constant Xmax- 

5. The function f'{w; 6) = {d/d9)f{w; 6) is absolutely integrable. 

6. For any given family, there is a positive probability of at least two failures. 

7. Defining 7r(s) = E[/(T,o < s) ^™!, we have 

/" \°{t] 
^^^ds<oo for all M G [0,r] and r = 1,2,3. (12) 

This assumption is needed in the analysis of the first-stage estimator. For r = 1, it 
parallels Assumption (5.4) of Keiding and Gill (1990), 

8. The matrix [{d / dj)\J , Ao{- , j))]\^=^o is invertible with probability going to 1 as 
n oo. 

9.2 Technical Preliminaries 

With 'ijj*{r,h) as in ([H]), we define il^mini^) = T^^T^o<r<m'ip*ir, h) and ipmaxW 
= msiXo<r<mi^*{r, h). In ([8]), the numerator and denominator are bounded above since 
Ui is assumed to have finite (m + 2)-th moment. Also, since Ui is nondegenerate, the 
numerator and denominator are strictly positive. Thus i'maxW is finite and 'ipminW is 
strictly positive. We present below two lemmas. The first lemma, which can be proved 
by elementary calculus, is taken from Zucker et al. (2006). The second lemma parallels 
Lemma 3 of Zucker et al. (2006). 

Lemma 1: The function tlj*{r,h) is decreasing in h. Hence for all 7 G ^ and all t, 
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ipi{'y,A,t) < V'max(O) i'iijy-^jt) > i^mini^^ ^i^)) ■ addition, there exist 5 > and 
/i > such that, for all h>h, ip^^„{h) > Bh^^ . 

Lemma 2: For any e > 0, we have sup^gj^ ,^] |Ao(s,7°) — Ao(s— ,7°)| — > as n oo. 
9.3 Consistency 

As indicated in Sec. 4, the consistency proof proceeds in several stages. 

Claim Al: Ao(t, 7) converges in probability to some function ^.^[t, 7) uniformly in t and 

7. The function A^(t,7) satisfies A*(t,7°) = A°(t). 

Remark: We give here an in pr. consistency result, rather than an a.s. result as in 
Zucker et al. (2006). The reason will be explained in the course of the proof. 
Proof: We can write Ao(t, 7) as 



Aoit,l)= / - ^ ^ . (13) 

Jo n-i Er=i nTro < s)Us-, 7, A) EJi Yrjis) exp(/3^Z,,) 

The proof here builds here on that of the corresponding Claim A in Zucker et al. The 

main point needing attention here is the fact that, because of the indicators /(Tjo < s), 

the denominator of f|T3l) tends to as s ^ 0. Special arguments are needed to deal with 

this "vanishing denominator" problem. 

Define, in parallel with Zucker et al. (2006), 



Sn(t,7,A) 

and 



n 



n-' Ell < 7, A) Er=i M^) exp(/3' Z,,) 



E[/(T,o < s)^,(s-,7,Ag) Er=i^^^-(^)^^P(/3 Z,,)] ^^o, 
E[/(r,o < s)^,(s-,7,A)Er=:i^ui^)exp(/3^Z,,)] ° 



By definition, Ao(t,7) satisfies the equation Ao(t, 7) = S„(t, 7, Ao(-, 7)). 

Remark: In Zucker et al. (2006), we had a result to the effect that ^^(t, 7,A) — >■ 
S(t,7, A) a.s. as n ^ 00, uniformly over t G [0,r], 7 G ^, and A in a certain set. We 
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could not obtain the corresponding result here; the argument of Aalen (1976) fails in the 
neighborhood of zero because of the vanishing denominator problem. This is why we give 
only an in pr. consistency result rather than an a.s. result. 
Again in parallel with Zucker et al. (2006), define 



^^^"'"^^ E[J(T,o < 5)^,(5-,7,A)E7=i^..(^)exp(rZ,,)] ^'^'^^ 
This function q~f{s, A) has the same properties as noted for the corresponding function in 
Zucker et al. These properties are not interfered with by the insertion of the indicator 
function /(Tjo < s). In particular, from Lemma 1 we have 



V'i(s-,7,A) ^^inihmax)' 

Pulling this bound outside of the expectation, we get the a bound on q~f{s, A) analogous to 
that in Zucker et al. Similarly, as in Zucker et al., the function q^{s,A) has the following 
Lipschitz-like property: \q^{s, Ai) — q^{s, A2)\ < K su.pQ^^^g\Ai{u) — A2{u)\. Accordingly, 
we find that the equation A{t) = S(t,7, A) has a unique solution, which we denote by 
Ao(t, 7). The claim then is that Ao{t, 7) converges in pr. (uniformly in t and 7) to Ao(t, 7). 
We now define, for any e > 0, the quantities 



S„(i,7,A,e) = ^ 



n-' Eti mo < sms-. 7, A) E7=i y^As) ^mi^' z^.) 

and 

. _ E[7(r,o < s)Us-, 7, Al) Y,^{s) exp(/3°^Z,,)] 
''^ Je E[7(r,o<«)V5.(«-,7,A)E;::i>^..(«)exp(/3%)] 

We next define Ao(i, 7, e) to be the solution of the equation Ao(t, 7, e) = 5„(t, 7, Ao(-, 7), e), 
starting from Ao(e, 7, e) = 0. We extend the definition of Ao(t, 7, e) by setting it equal to 
for t < e. Similarly, we define A^it, 7, e) to be the solution of the equation Ao(t, 7, e) — 
5(i, 7, Ao(-, 7), e), starting from Ao(e, 7, e) — 0, and extend the definition by setting 
AQ(t, 7, e) equal to for t < e. 
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For t e [e, r], the difference between AQ(i, 7) and AQ{t, 7, e) is as follows: AQ{t, 7) is the 
solution to Ao(t,7, e) = 7, Ao(-, 7), e), starting from Ao(e, 7,e) = Ao(e, 7), whereas 
AQ(i,7, e) is the solution to Ao(t,7, e) = H(i, 7, Ao(-, 7), e), starting from Ao(e, 7, e) — 0. 
Hence, by an induction argument similar to that in the proof of Hartman (1973, Thm. 
1.1), we find that 

|A*(i,7,e)-A*(t,7)| <e^Ao(e,7), 
where K is the Lipschitz constant for q^{s,A). We thus have 

sup |A*(t,7,e)-A*(t,7)|^0 as e ^ 0. (14) 
■yeg,te[o,T] 

Now, for any given e > 0, there is no vanishing denominator problem on the interval 
[e, r] . Hence, the argument in Zucker et al. (2006) goes through as is, and we get the 
following result: for any e > 0, 

sup |Ao(t, 7, e) — Ao(i,7, e)| — > a.s. as n — > 00. (15) 

•7Ge,te[e,r] 

In fact, in the supremum above, we can replace [e, r] by [0, r] , since by definition Ao(t, 7, e) = 
Ao(i,7, e) = for i < e. 

The above a.s. result immediately yields the corresponding in pr. result: 

sup |Ao(i, 7, e) — AQ(t, 7, e)| ^ in pr. as n ^ 00. (16) 

TGa,te[0,r] 

Our aim now is to show that 

sup |Ao(t, 7) — AQ(t,7)| ^0 in pr. as n ^ 00. (17) 

TG0,tG[O,r] 

That is, we want to show the following: for any p,d > 0, there exists n*(p, S) large enough 
such that 

Pr( sup \~Ao{tn)-Kit:7)\> p)<S 
-ye0,te[o,r] 

for all n > n*{p, 5). 
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Let p and 5 be given. By fll4p . we can find e > small enough such that 

sup |Ao(t,7,e) - AS(^,7)I < ^• 
Further, for this fixed e, the result (flGll implies that there exists n such that 

Pr( sup |Ao(t,7,e)-Ao(t,7,e)| > ^ I < ^ 

for all n>n. 
Now 

|Ao(t,7)-AS(t,7)l 
< |Ao(t, 7) - Ao(t, 7, e) I + |Ao(t, 7, e) - AS(t, 7, e) I + I AS(t, 7, e) - AS(t, 7) 



:i8) 



The developments just above imply that, for n > fi, the supremum over 'y E Q and 
t G [0, r] of the sum of the last two terms is bounded by |p with probability at least 
1 — |5. It remains to deal with the first term. 



Define 

Ci{s) = 

C2{S) = 

We can then write 



^ n nii 

-^/(Tio < s)7/;,(s-,7,A(-,7))^F,,(s)exp(/3%), 
^ i=i j=i 

-Y^I^T,, < s)^,(s-,7,A(-,7,e))^F,,(s)exp(/3'^Z,,). 



Ao(t, 7) - Ao(i, 7, e) = Ao(i A e, 7) + A{t, e). 



(19) 



where 



A{t,e) 



J the 



^ n rrii 
i=l 7 = 1 



We deal with the two terms on the right side of fll9l) in turn. In what follows, we let R 
denote a "generic" constant which may vary from one appearance to another, but does 
not depend on the unknown parameters or e. 
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Denote 



and recall the definition 7r(s) = E[7(rjo < s) Yl^li^iji^)]- ^^^^ recall 



iT^Umo < sM(s-,7,A(-,7))E7=:i^^.(5)exp(/3^Z,, 



It is clear that A(i,7) < RT{t,'y), where 



t 1 Y^n 



n(.) 



We can write 

T(t,7) 



/ n(s)-^ -V/(T,o<s) VF,,(s)exp(/3°X) 



Ao(s)(is 



+ 



where My is the martingale process corresponding to Nij: 

Mij{t) = Nij{t) - [ Xo{u) exp{(3°^Zij)Yij{u)ilJi{Y, K (20) 
Jo 

The first term is clearly bounded by RAQ{t). Thus, denoting the second term by 
M*{t), we have 

A{tn)<Rmt)+ sup \M*{u)\] (21) 

We next examine A{t, e). We can restrict to t > e, since A{t, e) = for t < e. Denote 
A{t) — Ao(i, 7) — Ao(t, 7, e). Bearing in mind the Lipschitz property of ■0, we find that 

\A{t,e)\<Rl\A{s-)\dT{s). 

Note that, for t > e, dA{t,e) = dA{t). Thus, a simple induction and some additional 
simple manipulations lead to the following, where we employ the symbol V to denote 
product integral and use the fact that A(e) = A(e, 7): 

\A(t,e)\ < \A(e)\Vl(l + RdT(s)) < |A(e)| exp(i?[T(t) - T(e)]) < |A(e, 7)! exp(i?T(T)) 
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In view of the analysis above of T{t), we get 



|A(t,6)|<|A(6,7)|exp(i?[A°(r)+ sup \M*{u)\]). 

ue[o,T] 



(22) 



Putting ([21]), and ([22]) together, we get 
|Ao(t,7)-Ao(t,7,e)l 



<i?i[A°(e)+ sup |il/r(M)|](l + exp(/?2[l+ sup \M*{u)\]). 

ue[0,T] M6[0,r] 



(23) 



M*(u) 



n(s 



-1 



for suitable absolute constants -Ri and R2. 

The last main step is to analyze the martingale process 

^ i=i j=i 

Our argument is pattered after the argument given by Keiding and Gill (1990, p. 595). 

By Lenglart's and Markov's inequalities, we have, for any positive k and r] and any 
ce [0,r], 

Pr(v^ sup \M*{t) \ > k) <r] + Pr(n(M*)(c) > tj/k^) < t] + —E[n{M*){c)]. 

u€[0,c] V 



Define J(s) = /(n(s) > 0). Then 



n(M*)(c) = / n(s) 
'0 



U{sr' - J2 <s)Y^ M^) exp(/3°^Z,, 
.'^ i=i j=i 



XQ{s)ds 



nJ{s) 



/o nU{s) ^ - Jo nU{s) + 1 



n + 1 



Ao(s)(is. 



As in Keiding and Gill, 



E 



n + 1 



nU{s) + 1 



< 



TT S 



Hence E[n(M*)(c)] < -R^i(c), where ^1 is as in Assumption 7. We thus get 



Pr( sup \M*{u)\ > K,n'^) < r] + RK^ri'^^i{c) 

mG[0,c] 



(24) 
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Now, the main quantities in the bound in are AQ(e) and sup„g[o,T] By decreas- 

ing e if necessary, we can make AQ(e) as small as we need. The behavior of sup^gjQ \M*{u) \ 
is characterized by i^^. We see that by decreasing e if necessary and choosing rj appro- 
priately, we can guarantee that the probability that the right side of (1231) is less than |p 
will be at least 1 — ^6 for all n sufficiently large. With this, we have taken care of the 
first term of f|T8|) . The desired convergence has thus been established. It is easy to see 
that A5(t,7°) =AS(t). 

Claim A2: Ao(t, 7) converges in probability to some function Ao(t, 7) uniformly in t and 
7. The function Ao(t,7) satisfies Ao(t, 7°) = ^o(t). 
Proof: We can write Ao(t,7) as 

Er=i ^^is- 7) J:T=i Ms) exp(/3^Z,,) ' 
In view of Claim Al above, up to a uniform error of op(l) we can replace all instances 
of Ao('U,7) in the definition of ipii^—,^) by Aq(m, 7). The desired result then can be 
obtained using the argument used to prove Claim A of Zucker et al. (2006). 
Claim A3: We have 

sup |Ao(s, 7°) — Ao(s— , 7°)| as n — > 00, 

se[0,T],76e 

sup |Ao(s, 7°) — Ao(s— , 7°)| as 00. 

sG[0,T],7e£; 

Proof: These results follow from Claims Al and A2 and the fact that AQ(t, 7) and Ao(t, 7) 
are continuous. 

Claim B: U(7, Ao(-, 7)) converges in probability uniformly in t and 7 to a limit u(7, Ao(-, 7)). 
Proof: As in Claim B of Zucker et al. (2006). 

Claim C: There exists a unique consistent (in pr.) root to U(7, Ao(-, 7)) = 0. 
Proof: By appeal to Foutz's (1977) theorem, as in Claim C of Zucker et al. (2006). 



Ao(t,7) 
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9.4 A workable representation of Ao{t) — A^it) 

In order to develop our asymptotic normality result, we need a workable representation 
of Ao(i) — A.Q(t). The first step is to develop a suitable representation of Ao(i) — A.^{t). 
Then, building on this, we develop our representation of Ao(t) — A.Q{t). 

9.4.1 Representation of Ao(t) - Ao(t) 

Our starting point is the following simple lemma. 
Lemma: Let Tln{t) and Sn{t) be stochastic processes, and let An{t, e) and -B„(i, e) be 
quantities that are bounded in probability uniformly in t and e. Define 

nn{t,e) = Unit) -An{t,e)nn{e), 
Sn{t,e) = Bn{t,e)[Sn{t)-Snie)]. 

Suppose that: 

p 

1. sup(g[g,-] ^yn\^Zn{t, e) — Sn{t, e)| ^ as n ^ oo for any fixed e > 0. 

2. lim4olinisup„_^(^Pr(sup^g[o,e] V^l^n(^)| > 5) = for all S > 0. 

3. lim£j^oliKisup„^3oPr(supj£[o,e] V^\^n{t)\ > 5) = for all d > 0. 

4. Bn{t, e) — > Bn{t, 0) uniformly in t as e — > with probability converging to one as 
n — > oo. 

Then sup,e[o,r] V^l^n(t) - B^{t, 0)5„(t, 0)| ^ 0. 

We apply this lemma with TZn{t) = -\/n[Ao(^) — AQ(i)]. We have to check the four 
conditions enumerated in the lemma. 

Cuuditioii 1 
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Arguments along the lines of Zucker et al. (2006) yield the result of Condition 1, with 

2n rrii 



ft ~( \ 1 "H 



i=i j=i 



(25) 



where 



pit,e)= n 

s&[e,t] 



2n nii 



Here 



l + n-^^^aj(s,t)(i%(s) + n-^0*(s)/(rio < s)S, 

i=l j=0 



^ 2n 



1=1 



with Ri.(s) = E^Ji eMf3°^'Zij): 



k=l 



1=1 



Ri.{u)rju{u) exp{l3°'^Z 



iO 



k=l 1=1 



and for j > 1 



(26) 



%{s, t) = n / 2^ 2^ /(Tfco < M)o?iVfci(M). 

-/s xyy^^ ^q)s k=i 1=1 



In the above, r]ii{s) is defined as 

^ 03i(7°,Ag,g) _ f Ml\Ks) 

In Sec. 9.3.2 below, wc present in detail a similar argument for Ao(t) — ^^{t). 

Appealing to Assumption 7 and using arguments similar to those used in the consis- 
tency proof, we find that the Q quantities defined above converge in probability uniformly 
in s and so that p(i, e) converges in probability to a deterministic limit uniformly in t 
and e. 

Condition 2, 3, and 4 
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In regard to Condition 2, we have 



Ao(t,7)-AoW=Ai(t) + A2(t), 



where 



A2(t) 



n-'Eliimo<s)E7:idM,,{s) 



where 



r(s,7) 



n-' TJU mo < s)Us-, 7, A) Er=i ^^As) exp(/3^Z,,) 



n-' Er=i ^7^.0 < 7°, Ag) Er=i exp(/3°^Z 

n-' Eti m^ < 7, A) ET=i M^) exp(/3^Z,, 



and Mij{t) is defined as in (l20l) . We will deal with Ai(i(:) and A2(t) in turn, starting with 
A2(t). In the development below, R denotes a "generic" absolute constant. 
The quadratic variation process of A2(t) is given by 



(A2)(t) 



n-'J:tlnT^o < s)^.(g-,7°,Ag)Er=i>^ij(g)exp(/3°%: 
[n-' Ell HTio < s)^.(s-,7,A) E7=l'^^A')^Mf3^^^JW 



Ao(s)(is. 



By arguments similar to those used in connection with M*{t) in the proof of Claim Al, 
we find that E[n(A2)(t)] < RC,i{t). An application of Lenglart's inequality then gives 

PT{y/n sup |A2(t)| > k) <r] + RK^7]-^ii{t) Wr] > 0. 

iG[0,e] 

Assumption 7 implies that ^i(e) | as e | 0, and this takes care of A2(t). 

We now turn to Ai(t). As before, denote J(s) = /(n(s) > 0). We can write 



with 



and 



A,a{t)= [\T{s,-f)-l]J{s)X°{s)ds 
Jo 



Ai6(t) = / [As) - l]X°is)ds. 
'o 
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The term A.ib{t) can be shown to be uniformly Op{n 2) by the argument in the middle 
of page 595 in Keiding and Gill (1990) As for Aia{t), we have 

Aiait) < Rt\Ao{t,-f)-A°{t)\ < Rt\A,{t)\+Rt\A2it)\ < Rt\A,ait)\+Rt\Auit)\+Rt\A2it)\. 
Thus, for t small, 

Rf 

and the terms on the right hand side have already been taken care of. 

The proof of Condition 3 is similar to that given above for A2{t). Condition 4 follows 
easily from the uniform convergence of the quantities. 

9.4.2 Representation of Ao{t) - Ao{t) 
Let 

2n 

y{s, {Ao,Ao}) = -J2MY,Ao,s)RiXs) 
1=1 

and 

y{s,A) = -y2^,{Y,A,s)Rds) 

1=1 

so that in Aq, s) we take Ao(Tjo) if Tjo > s and Ao(Tjj) if Tj^ < s, j > 0. By Claim 

A3, we have that also sup^gjo^^j |Ao(s, 7°) — Ao(s— , 7°)| converges to zero. Thus, we obtain 
the following approximation, uniformly over t e [0, r]: 

Ao(t,7°)-AoW«- / {y(s,K)}-'^J2dM,,{s) 

^ -^0 i=l j=l 

+- / \{y{sd~^o,Ao})}-'-{y{s,Al)}-']j2T.^N,,{s). 
^ -^0 ^ i=l j=l 

Now let 



*M -1 



Xis,r)^{yis,Al + rA*)} 

with A* = Aq — Aq or Aq — Aq, according to the estimator being used. Define X and X 
as the first and second derivative of X with respect to r, respectively. Then, by a first 
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order Taylor expansion of X{s, r) we get 

i=l j=l 



/ E Toni E ^(^« ^ ^) exp(/3°^Z,,){Ao(r«) - A°M} E ^^^^■(^) 
/ E ;v .T°n2^ (^^° ^ ^) exp(/3°^Z,o){Ao(T.o) - Ao(r.o)} E E ^^iAs) 
' ' E < ^) exp(^°^Z,o){Ao(T,o) - A°(T,o)} E E ^^^^•(^)- 



— n 



The justification for ignoring the remainder term in the Taylor expansion is as in the 
parallel argument in Zucker et al. (2006). 

The second, third and fifth terms of the above equation can be written, by inter- 
changing the order of integration, as 

2n mi 



-n 



„j ^,1 1,1,1 

/ {Ao(5)-A°(5)}^^T,,(5,t)rfiV,,(5) 

•^0 i=l j=0 



where 

T f .^ -1 Ri.{u)riii{u) exp(/3°^Z^o) ^ ^ 

= " / was)? £ S """^"^ 

and for j > 1 



Tjj(s,t) = n 



-1 / Ri.{u)T]u{u) exp(/3° Z^j) . . 



The fourth term can be written, by plugging in the representation for Aq — Aq, as 



Jo y[S, Aq) .^^ .^-^ 
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where 



A{s, t) = r X; ?#f# exp(/3°^Z,o) [ F {p{v)}-'dNM] ^ f; dAr^,(.) 
and -/Vfeo(^) ~ -^(^feo ^ Given all the above, we get 

•^0 i=l 7 = 1 



2n mi 



A{s,t)p{s-] 



J]J]/(r,o<.)dM,,.(.) 



Jo 



-n-^ / {ko{s)-h.l{s)}Y,Y.^i3MdN,^{s) 

-^0 i=i ,=0 



2n mi 



By solving the above approximation recursively, for the relatives' failure times, we get 



Ao(t,7°)-AoW 



p{s- 



2n mi 



np{t) „ 



2rt rrii 



p{t-)dN{t) r 



2n irii 



1 yy^^ -'^oJ 



where A^(t) = E?:iEr=\^^.W, 



1=1 J 



2n nii 



B(s, t) = n^' I A(s, u)p{u-) V V Tijiu, t-)dNij(u) 
'o 



i=l j=Q 



and 



p(^) = n 



2n m: 



i=l j=0 



9.5 Asymptotic normality 



To show that 7 is asymptotically normally distributed, we write 



0= U(7,Ao(-,7)) 
= U(7°, Ao) + [U(7°, Ao(-, 7°)) - U(7°, A°)] 
+ [U(7,Ao(-,7))-U(7°,Ao(-,7°))]- 
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We examine in turn each of the terms on the right-hand side of the above equation. 
Step I 

We can write U(7°, Aq) as 

(n 2n \ 

i=l i=l / 

Here i — 1, . . . ,n are iid mean-zero random (p -|- l)-vectors stemming from the hkeh- 

(2) 

hood of the proband data, while i = 1, . . . ,2n are iid mean-zero random (p+l)-vectors 
stemming from the hkehhood of the relatives' data. It follows immediately from the clas- 
sical central limit theorem that n~^/^U(7°, Aq) is asymptotically mean-zero multivariate 
normal. 
Step II 

Let Ur — Ur{'y°, Aq), r — 1, . . . ,p, and Up+i — C/p+i(7°, Aq) (in this segment of the 
proof, when we write (7°,Ao) the intent is to signify (7°, Ao(-, 7°)). First order Taylor 
expansion of Ur about Aq, r = l,...,p-|-l, gives 

n'/'{Urh°,Ao)-Ur{Y,K)} 

2n mi 

= n-'/'J2IlQ^^r{l°,KTij){MTij,Y) - KiTij)} + Op(l), (27) 

i=l j=0 

where 

OU f'Y° A°) 

QijriY, Ao, Tij) = g^o^^,,^ i = l,...,2n j = 0, . . . , r = 1, . . . ,p + 1. 
Now let N^j{t) = I{Tij <t) i = l,...,2n J = 0,...,mi. Then 

^ 2?T. nT'i iff 

Ur{Y,h)-Ur{Y,K) = -J2Yl / QiJr{Y,A°,s){ko{s,Y) - Kis)}dN:^{s). (28) 

^ i=i j=o Jo 

Define 
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and 

y{u,K)nj-ijr'Jo Pis) 

Also, let ar"^''(-u), Q;r^'*(-u) and al?\u) denote the corresponding limiting values of Ai^\u), 
Ar (u) and Ar (u) as n goes to infinity. Then, after plugging into ( !28l) the representation 
in Sec. 9.3.2 for y/n[Ao{s,^°) — Aq(s)] and replacing the ^'s with their limiting values, 
we obtain 

t/,(7°,Ao)-f/.(7°,A°) 

fc=l /=i 

This gives a representation of f/r(7°, Aq) — ?7r(7°, Aq) r = l,...,p+las the average of 
independent mean zero iid random variables. Hence, asymptotic normality follows from 
the classical central limit theorem. 
Step III 

First order Taylor expansion of U(7, Ao(-,7)) about 7° = {(3°'^, 6°)'^ gives 
U(7, Ao(-, 7)) = U(7°, Ao(-, 7°)) + D(7°)(7 - 7°)^ + o,{l), 

where 

DUl)=dUi{j,Aoi;-f))/d^, 

for /, s = 1, . . . ,p + 1. 

Combining the results of Steps I-III above we get that 71,^/^(7 — 7°) is asymptoti- 
cally zero-mean normally distributed with a covariance matrix that can be consistently 
estimated by a sandwich-type estimator. 
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Table 1: Simulation results: 500 control proband matched with 500 case probands; One 
relative for each proband; (3 — 0.693, Ao(i) — t , 9 — 2.0, 500 samples. 





Proposed Method 




Hsu et al. 


Shih and Chatter jjee 






Empirical 




Empirical 




Empirical 




mean 


Standard Error 


mean 


Standard Error 


mean 


Standard Error 


^ 


0.706 


0.197 


0.697 


0.201 


0.698 


0.182 


9 


2.003 


0.312 


1.986. 


0.302 


1.992 


0.303 


Ao(0.2) 


0.201 


0.034 


0.204 


0.030 


0.202 


0.029 


Ao(0.4) 


0.402 


0.063 


0.407 


0.058 


0.403 


0.054 


Ao(0.6) 


0.603 


0.095 


0.612 


0.090 


0.605 


0.084 


Ao(0.8) 


0.809 


0.136 


0.820 


0.131 


0.811 


0.122 
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Table 2: Simulation results for the proposed estimators: 500 control proband matched with 
500 case probands; One relative for each proband; Ao(t) = t; 90% censoring rate; 2000 
samples. 
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Table 3: Simulation results for the proposed estimators: 500 control proband matched with 
500 case probands; One relative for each proband; Ao(t) = t; 60% censoring rate; 2000 
samples. 
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Table 4: Simulation results for the proposed estimators: 500 control proband matched with 
500 case probands; One relative for each proband; Ao(t) = t; 30% censoring rate; 2000 
samples. 
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Table 5: Simulation results of left-restricted data: 500 control proband matched with 500 
case probands; One relative for each proband; sq — Q.l, (5 — 0.693, h.Q{t) — t , 9 — 2.0; 
500 samples. 
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Table 6: Analysis of a case-control family study of breast cancer. 
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